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ABSTRACT 

The flow pattern and solidification 
phenomena in die casting process have been 
investigated in the first phase study. The flow 
pattern in filling process is predicted by using a 
VOF method. A good agreement with experimental 
observation is obtained for filling the water into a 
die cavity with different gate geometry and with an 
obstacle in the cavity. An enthalpy method has 
been applied to solve the solidification problem. 
By treating the latent heat implicitly into the 
enthalpy instead of explicitly into the source term, 
the CPU time can be reduced at least 20 times. The 
effect of material properties on solidification fronts 
is tested. It concludes that the dependence of 
properties on temperature is significant. The 
influence of the natural convection over the 
diffusion has also been studied. The result show's 
that the liquid metal solidification phenomena is 
diffusion dominant, and the natural convection can 
affect the shape of the interface. In the second 
phase study, the filling and solidification processes 
will be considered simultaneously. 

INTRODUCTION 

Due to its importance in manufacturing, 
analysis of casting process for metals has received 
considerable attention in recent years. For 
improving the quality of manufactural parts with 
minimal experimentation, it is important to 
understand how the transport phenomena could 
affect the microstructure. Numerical analysis and 
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modeling efforts provide a powerful means to 
achieve the goal. 

For simplifying the problem, die casting 
process can be basically divided into the filling 
stage and solidification stage. It can be assumed 
that no phase change in the filling stage and no 
forced flow when phase change started, because 
the filling is much faster than solidification. The 
study of both of these two phenomena is presented 
in this paper, and the combined process will be 
investigated in the next step of the study. 

The fluid flow in filling stage involves a free 
surface between two different fluids (gas and 
liquid). The volume of fluid (VOF) method 
developed in SOLA-VOF 1 * is used here and the 
numerical implementation details can be found in 
Ref [2]. The surface tension force is calculated by 
Continuum Surface Force (CSF) model 3 , in which 
the surface tension is treated as a body force. 

The solidification is handled by an enthalpy 
model, which is initially developed by Shamsunder 
and Sparrow 4 . In this model, the enthalpy is used 
as a dependent variable. The liquid phase, solid 
phase and the interface all satisfy a unique 
governing equation. This can be solved like a 
single phase problem, and no interface tracking 
necessary. After the enthalpy field is solved at 
every time step, the temperature and the interface 
position can be obtained from the enthalpy- 
temperature diagram. Voller and coworkers 5 * 7 
improved this model by applying a porosity 
concept to identify the phase change region, which 
forms a so-called enthalpy-porosity model. In this 
model, the latent heat content, AH, varies among 0 
(cell all solid) and L (cell all liquid), w'here L is the 
latent heat of the phase change. The porosity X is 
defined as the ratio of AH over L, then X decrease 
from 1 to 0 when phase changes from liquid to 
solid. So, by prescribing a “Darcy” source term in 
the momentum equation, w'hich vanishes when X is 
0 and increases rapidly when X approaches to 1, 
the velocity will be forced to 0 when total 
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solidification is achieved. All the previous 
calculations in Ref. [5-7] put the latent heat term 
explicitly into the source. In turn, the iteration is 
necessary in every time step for updating the latent 
heat content, so that the temperature in the cell that 
is undertaking phase change will keep its melting 
temperature. As the experience of the authors and 
as pointed in Ref [6], the average iteration times is 
about 30-40 per time step. This is very time 
consuming and expensive, which will make the 3- 
D simulation unpracticable. In this paper, the latent 
heat will be included in the enthalpy and no 
iteration is needed, while the CPU time will 
decrease at least 20 times with the only extra effort 
of getting the temperature from enthalpy. 


GOVERNING EQUATIONS 

The die casting process can be assumed as a 
Newtonian fluid laminar flow. No turbulence is 
considered because the duration time of this 
transient process is very short. The governing 
equations in 2-D Cartesian coordinates are as 
follows. 

Continuity equation: 
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Momentum equation in y-direction: 
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be include in the momentum equations. It will be 
discussed in the CSF model. The S f term is gravity 
force -pg and the “Darcy” source terms S x and S y 
are zero. The energy equation does not need to be 
solved for this isothermal process. 

In the solidification stage, all the properties, 
including density, viscosity and thermal 
conductivity, are constants when using the 
incompressible flow approach. The natural 
convection is generated by the buoyancy force 
using Boussinesq approximation. 

S f = pg/3(T-T ref ) (5) 

Where, p is thermal expansion coefficient 
and T rcf is the reference temperature. Because the 
temperature range is very large in solid phase 
(from room temperature to melting temperature), 
the dependence of properties on temperature can 
not be ignored. In this case, the density, viscosity, 
thermal conductivity and heat capacity will be 
updated on every time step according to the 
temperature. The S f term will just be a gravity 
force -pg. The “Darcy” source terms S x , S y and the 
energy source term S h will be given later in the 
enthalpy-porosity model. 


VOF METHOD and CSF MODEL 


In the filling stage of casting process, the 
position and shape of the free surface between the 
liquid and gas phase can be represented by the 
VOF values. The fractional volume in a typical 
control volume cell is defmed as: 




V +V, 
g I 


( 6 ) 


Where, V represents volume occupied by gas (V g ) 
and liquid (V L ) within the control volume 
considered. The function F obeys the volume flux 
conservation equation: 



du F 

^7 


= 0 


The average density is then defined as 


(7) 


In filling stage, the free surface exists and 
the surface tension force terms F sx and F sy should 


p= Pl F+p g {l-F) (8) 
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The surface tension is an inherent 
characteristic of material interface, because fluid 
molecules at or near the surface experience uneven 
molecular forces of attraction. CSF model 
interprets the surface tension force as a continuous, 
three-dimensional effect across an interface, and 
the interface reconstruction is no longer needed. 


F s = okVF (9) 


Where cr is the fluid surface tension coefficient (in 
units of force per unit length), F is the VOF value 
and k is the local surface curvature. 


k = 


-V 


r vf'' 


ivn 


( 10 ) 


ENTHALPY-POROSITY MODEL 


Voller 5 ' 7 improved the enthalpy model by 
introducing the porosity concept, which forms the 
enthalpy-porosity model. In this model, an unique 
energy equation is suitable for the whole 
calculation domain and no liquid-solid interface 
need to be accounted for. The zero velocity 
condition in the vicinity of the interface is modeled 
by a “Darcy” source term in momentum equation. 

Voller substitutes the H=cT+AH into the 
energy equation and moves the time variation and 
convection term of AH to the right hand side, 
which forms (relate to Equation (4)) 


H = cT (11) 

dpt^H d 

S,=— + - (pu .AH) (12) 

h a a* J 

J 

In the approach of this paper, however, AH 
is implicitly included in the enthalpy so that 


H = cT + AH 

(13) 

S h = ° 

(14) 


After the enthalpy field is obtained, the 
temperature and the porosity can be solved from 
enthalpy-temperature diagram: 


IfT<T m -s, 

then H = c T, AH = 0 
s * 

IfT - s <T <T +£, 

J m m ' 

then H ~ c T + A//, 0<A H<L (15) 

then H = cT+ L + Cj{r-^T m + e )J, 

AH = L 


The subscripts s, 1 denote solid phase and 
liquid phase respectively. If the material is not pure 
(e.g. alloy), the phase change takes place over a 
temperature range, call mushy region (T ra -s to 
T m +e), where T m refers to the center of the mushy 
region. For isothermal phase change problem (pure 
material), T m is the melting temperature and e 
equal to zero. The problem considered in this study 
is isothermal phase change. 

The “Darcy” source terms in momentum 
equations are constructed as: 

S x =Au (16) 


Where 


S y= Av 

D(l-A ) 2 


(17) 


(18) 


The porosity X is: 


A = 


A H 
L 


(19) 


The constant D is very large that can force the 
velocity effectively to zero when the cell is total 
solidified. Constant d is a small number for 
avoiding divided by zero. In this study, D and d are 
chosen as D=1.6xl0 5 and d=0.1 according to Ref 
[5]. The parameter A is formulated in such a way 
that the “Darcy” source term has no any influence 
in the fully liquid cells. While it gradually increase, 
dominate and finally override all other terms in the 
momentum equation when solid portion increases 
in that cell. This smooth transition is superior than 
the step changes that may cause divergence. 


NUMERICAL METHOD 
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The flow solver used in this study is the 
FDNS code (Finite Difference Navier Stokes 
Solver). FDNS code can be used to solve 2-D 
planner, 2-D axisymmetric or 3-D forms of Navier- 
Stokes equations and other scale transport 
equations using primitive variables and curvilinear 
coordinate with non-staggered mesh systems. A 
pressure based predictor/multi-corrector solution 
procedure is employed in the FDNS code to 
enhance velocity-pressure coupling and mass- 
conserved solution at the end of each time step. 
This pressure based method is suitable for all speed 
flow computations including supersonic and 
subsonic flow, compressible and incompressible 
flow, laminar and turbulent flow. The detailed 
description of FDNS code can be found in Ref. [8, 

9]. 

RESULTS and DISCUSSIONS 

The water experiments 10 ’ 11 during the filling 
stage are calculated here for flow pattern 
comparison. The dimension of the die is 
152x102x3 mm with the gates and obstacle having 
the same thickness as the die (3mm). A uniform 
61x61 grid mesh is used in all the calculations. The 
gravity force is considered downwards. As pointed 
in Ref. [7] that due to gate velocity in the 
experiment may vary with time and the calculation 
assume constant inlet velocity, the comparison of 
the prediction and the experiment is made at the 
same percentage filling of cavity. 

Figure 1 is for one gate case, where the gate 
is located on 3 mm from the bottom of the right 
wall. The gate width is 6 mm, and the gate velocity 
is 10.7 m/s which corresponding to the Reynolds 
number of 6.4x1 0 4 . The void in the left bottom 
comer is caused by the recirculation and 
reattachment of the fluid on the bottom wall. A 
very fme grid mesh is needed to capture this void 
during the whole process. 

Figure 2 is for the case of two inflow gates, 
where one gate is at the bottom right wall and 
another is at the bottom left wall. The gate width is 
12 mm, and the gate velocity is 15.8 m/s which 
corresponding to the Reynolds number of 1.9xl0 5 . 
The two incoming flows impinge on and go 
upward, and then travel along the walls. The 
numerical predictions are verified by the 
experiments. 

Figure 3 is the case with an obstacle inside 
the cavity, where the gate is at center of the bottom 


wall. The gate width is 12 mm, and the gate 
velocity is 18.3 m/s which corresponding to the 
Reynolds number of 2.2x1 0 5 . The obstacle is 
50x25 mm located 25 mm above the gate. The 
figure illustrates that two recirculation zones are 
generated and preserved near the bottom left and 
right comers. Other two separation zones are 
present at the left and right walls of obstacle near 
the end of filling. The calculated flow patterns 
agree with the experimental observations very 
well. 


Table L Thermophysical properties of aluminum 


P (kg/m J ) 

2.38xlO J 

p (kg/m s) 

2.9x10° 

P (1/K) 

1.16x10“* 

c (J/kg K) 

1.08xl0 i 

k (W/m K) 

1.03x10* 

L (J/kg) 

4xKF 


Figures 4-7 represent the results of 
solidification of liquid aluminum. The liquid 
aluminum is initially 981.7 K (50 K superheat) 
inside a 150x150 mm cavity. At the beginning of 
the solidification, the left wall is lowered to 
291.7 K (640 K below melting point) and other 
three walls are adiabatic. Figure 4 shows the 
moving interfaces and liquid flow patterns for the 
constant properties approach, where the properties 
are given in Table 1. Case (a) considered the 
diffusion only and case (b) involved the natural 
convection also. Figure 5 is for the variable 
properties approach, where the density, viscosity, 
thermal conductivity and heat capacity are updated 
according to the temperature at every time step. 
Same as in Figure 4, case (a) considered the 
diffusion only and case (b) involved the natural 
convection also. Comparing case (a) with case (b) 
in Figures 4 and 5, one can see that the conduction 
dominates the solidification process, and the 
natural convection acts on it by effecting the shape 
of the interface. Because the hot fluid goes upward 
and cold fluid goes downward, the solid front goes 
faster at the bottom than it does at the top. The 
comparison of Figure 4 and Figure 5 shows that the 
dependence of the material properties on 
temperature cannot be ignored. For more clear 
demonstrating this influence, the positions of solid 
front at the bottom wall for case (b) are drawn in 
Figure 6. The difference can reach 50% for this 
particular case. The main reason for such a large 
effect is that the temperature changes 640 K inside 
the solid phase, so that the consideration of the 
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thermal conductivity and heat capacity as functions 
of temperature will largely enhance the heat 
conduction. 

Except accounting for the property effects, 
the main contribution of this paper is that one 
algorithm is developed to solve the energy 
equation by including the latent heat inside the 
enthalpy instead of in source as did in Ref. [5-7]. 
Because the iteration is not needed in each time 
step, this method saves CPU time by more than 20 
times. As shown in Figure 7, the correctness of this 
method has been justified by comparing the 
interface positions carried out by this method with 
the theoritical predictions 12 for the case of purely 
conduction with constant properties. When 
considering the convection, the phase front has 
been compared with the calculation in Ref. [7], as 
shown in Figure 6. These two methods got same 
solid front position at 28 second after the 
beginning of solidification. 

CONCLUSIONS 

In the liquid filling analysis, flow patterns of 
different gate position and with or without inside 
obstacle have been compared to the experimental 
observations. The good agreement shows that the 
numerical predictions can provide useful 
information of the filling patterns and the order of 
different areas in cavity to be filled. The accurate 
predictions of the filling stage establish a basic 
foundation for further resolving the solidification 
and finally the die casting process. 

In the investigation of the solidification, the effect 
of the variation of thermophysical properties as 
functions of temperature has been tested. This 
effect is very large due to the temperature range is 
very big for liquid metal solidification applications. 
It greatly reduces CPU time if solving energy 
equation by including the latent heat into the 
enthalpy instead of put it in the source term. This 
makes the simulation of 3-D real applications 
possible. 

REFERENCE 

1. Nicholes, B. D., Hirt, C. W., and Hotchkiss, R. 
S.: “SOLA-VOF: A Solution Algorithm for 
Transient Fluid Flow with Incompressible 
Flows w ith Free Surface,” Los Alamos National 
Lab., LA-8355, Aug., 1980. 


2. Chen, Y. S., Shang, H. M., Liaw, P., Chen, C. 
P., and Wang, T. S.: “ A Unified Two-Phase 
Numerical Method for General Gas-Liquid 
Flow Applications,” Gas-Liquid Flows, FED- 
Vol. 25, pp. 99-106, Symposium of the Joint 
ASME/JSME Fluids Engineering Conference, 
August 13-18, 1995, Hilton Heat Island, SC. 

3. Backbill, J. U., Kothe, D. B., and Zemach, C.: 
“A Continuum Method for Modeling Surface 
Tension,” Journal of Computational Physics, 
Vol. 100, pp. 335-354, 1992. 

4. Shamsundar, N., and Sparrow, E. M.: “Analysis 
of Multidimensional Conduction Phase Change 
Via the Enthalpy Model,” Trans, of the ASME, 
Journal of Heat Transfer, pp. 333-340, Aug. 
1975. 

5. Voller, V. R., and Prakash: “A fixed grid 

numerical modeling methodology for 

convection-diffusion mushy region phase- 
change problems,” International Journal of Heat 
and Mass Transfer, Vol. 30, No. 8, pp. 1709- 
1719, 1987. 

6. Brent, A. D., Voller, V. R., and Reid, K. J.: 

“Enthalpy-Porosity Technique for Modeling 
Convection-Diffusion Phase Change: 

Application to the Melting of A Pure Metal,” 
Numerical Heat Transfer, Vol. 13, pp. 297-318, 
1988. 

7. Minaie, B., Stelson, K. A., and Voller, V. R.: 
“Analysis of Flow Patterns and Solidification 
Phenomena in the Die Casting Process,” Trans, 
of the ASME, Journal of Heat Transfer, Vol. 
1 13, pp. 296-302, July, 1991. 

8. Chen, Y. S.: “FDNS: A General Purpose CFD 
Code, User’s Guide,” Engineering Sciences, 
Inc., Huntsville, Alabama. 

9. Chen, Y. S., Liaw, P., and Shang, H. M.: 
“Numerical Analysis of Complex Internal and 
External Viscous Flows with A Second-Order 
Pressure-Based Method,” AIAA Paper 93- 
2966, 1993. 

10. Smith, W. E., and Wallace, J. F.: “Gating of 
Die Castings,” Trans, of AFS, Vol. 71, pp. 325- 
348, 1963. 

11. Stuhrke, W. F., and Wallace, J. F.: “Gating of 
Die Castings,” Trans, of AFS, Vol. 72, pp.374- 
407, 1964. 

12. Gebhart, Benjamin: Heat Conduction and Mass 
Diffusion, McGraw-Hill Series in Mechanical 
Engineering, McGraw-Hill, Inc., 1993. 


5 

American Institute of Aeronautics and Astronaustics 




Figure 1 : Comparisons between photographs of experiments and calculated flow patterns with vector plots of velocity field 

for the casting with the gate located near the bottom of the right wall 







Figure 2: Comparisons between photographs of experiments and calculated flow patterns with vector plots of velocity field 
for the casting with the gates located at the bottom of the left and right walls 




Figure 3: Comparisons between photographs of experiments and calculated flow patterns with vector plots of velocity held 
for the casting with an internal obstacle and the gate located at die center of the bottom wall 






15 s t = 28 s 



(b). heat diffusion and natural convection 





28 s 



(b). heat diffusion and natural convection 









Time (s) 

Figure 6: Comparison of the solid-liquid interface positions for liquid aluminum solidification 
process, both conduction and natural convection are considered 



Time (s) 

Figure 7: Comparison of the solid-liquid interface positions for liquid aluminum solidification 
process, constant properties and conduction only 




